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The  purpose  of  this  project  was  to  develop  a  numerical  method  capable 
of  calculating  the  laminar  flow  about  a  two  dimensional  accelerating  body 
using  the  incompressible  Navier-Stokes  equations.  This  problem  arises  in  the 
calculation  of  the  dynamic  stall  of  an  airfoil  and  in  the  calculation  of  pitching 
and  heaving  stability  derivatives  for  an  isolated  airfoil.  It  is  also  a  first  step 
toward  solving  the  three  dimensional  flow  about  a  complete  wing  undergoing 
arbitrary  motion  at  high  angles  of  attack,  including  departure,  spin,  and 
post-stall  maneuvers. 

The  contravariant  form  of  the  momentum  and  continuity  equations  are 
derived  and  a  finite  difference  approximation  developed.  The  contravariant 
velocity  formulation  is  based  on  the  use  of  both  inertial  and  body-axis 
velocities  to  achieve  a  form  for  the  momentum  equations  that  has  the 
potential  for  improved  numerical  characteristics  for  rotating  coordinate  frames. 

Numerical  calulations  to  date  have  not  been  satisfactory.  Test  cases 
presented  are  Stoke’s  first  problem,  a  circular  cylinder  translating  at  a 
Reynold's  number  of  20,  and  a  rotating  stationary  circular  cylinder. 
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The  aerodynamics  of  aircraft  in  steady  flight  at  low  to  moderate  angles 
of  attack  are  veil  understood.  However,  the  viscous  flow  about  wings  and 
bodies  undergoing  arbitrary  unsteady  motion  is  not.  This  is  especially  true  for 
wings  at  high  angles  of  attack.  There  is  a  serious  need  to  be  able  to  predict 
the  performance  and  flying  qualities  of  aircraft  maneuvering  at  high  angles  of 
attack,  and  to  be  able  to  predict  departure  and  spin  behaviour,  and  to  exploit 
post  stall  maneuverability.  These  predictions  require  an  accurate  calculation  of 
the  nonlinear  forces  and  moments  that  are  associated  with  the  separated  flow 
under  such  conditions.  A  basic  requirement  is  to  be  able  to  calculate  the  flow 
about  an  isolated  wing,  and  the  analysis  of  a  simple  airfoil  is  the  first  step. 

The  aerodynamics  of  an  airfoil  which  is  pitching  are  influenced  by  the 
fact  that  the  local  velocity  of  different  points  on  the  airfoil  is  not  the  same  as 
the  velocity  of  the  reference  center.  This  has  implications  for  both  the  flow 
outside  the  boundary  layer  and  the  resulting  pressure  distribution,  and  for  the 
velocity  profile  within  the  boundary  layer  itself.  Vortices  that  are  shed  into 
the  flow  are  swept  downstream  and  cause  unsteady  effects  if  their  production 
is  not  continuous  and  constant.  In  addition,  large  separated  regions  exist 
which  may  vary  with  time.  The  convection  of  vortices  and  the  unsteady 
separated  regions  cause  time  lags  in  the  dynamic  forces  on  the  airfoil,  even 
for  incompressible  flow,  due  to  the  time  required  for  vortices  to  be  convected 
away  from  the  vicinity  of  the  body.  The  flow  at  the  body  surface  must  follow 
the  motion  of  the  body  in  order  to  satisfy  the  uo  slip  bounday  condition. 
Thus,  for  even  a  steady  external  flow  that  is  similar  to  the  unsteady  flow,  the 
boundary  layer  will  be  significantly  different,  because  the  flow  must  transition 
from  the  external  conditions  to  the  new  unsteady  and  spacially  varying 
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conditions  at  the  wall 

The  complexity  of  these  unsteady  phenomena  provides  a  strong 
incentive  to  use  numerical  methods  to  predict  these  flovs.  A  full  simulation 
of  the  three  dimensional  unsteady,  viscous,  and  compressible  effects  is  beyond 
the  capability  of  the  currently  available  methods  and  computing  resources. 
Therefore,  the  scope  of  the  problem  has  been  reduced  for  this  project  by 
neglecting  the  effects  of  compressiblity  and  turbulence,  and  by  restricting  the 
analysis  to  two  dimensional  flow. 

Several  investigators  have  used  numerical  methods  to  calculate  the 
unsteady  flov  about  bodies.  Hodge  (Ref  1)  considered  an  airfoil  in  steady 
motion  at  high  angles  of  attack.  The  flov  under  these  conditions  (ten  degrees 
angle  of  attack,  a  Reynolds  number  of  200000,  and  laminar  flov)  vas 
unsteady,  vith  separation  bubbles  forming  and  collapsing  on  the  upper  surface. 
The  numerical  technique  used  vas  vey  similar  to  the  present  method. 
Primitive  variables  vere  used.  The  finite  difference  approximation  used  a  first 
order  backvard  difference  in  time,  three  point  one  sided  upvind  differences 
vere  used  for  first  derivatives,  and  central  differences  vere  used  for  second 
order  derivatives  and  pressure  gradients.  The  equations  vere  solved  through 
successive-over-reiaxation  (SOR).  The  Poisson  equation  vas  used  for  pressure 
calculations.  The  rate  of  convergence  of  the  calculations  vas  driven  by  the 
convergence  of  the  pressure  on  the  vail  at  the  sharp  trailing  edge.  This 
experience  provided  the  motivation  for  the  present  cell  centered  pressures. 

Hegna  (Ref  2)  extended  Hodge’s  vork  by  adding  a  turbulence  model  and 
considered  a  pitching  airfoil.  Fictitious  body  forces  vere  used  to  account  for 
the  accelerations  due  to  the  rotating  reference  frame.  Constant,  linearly 
varying,  and  sinusoidal  pitching  motions  vere  considered. 

Mehta  (Ref  3),  used  a  stream  function-vorticity  formulation  to  solve  for 
the  flov  about  an  airfoil  oscillating  in  pitch.  The  vorticity  equation  vas  based 
on  the  body-axis  equations  vith  fictitious  body  forces. 


More  recently,  Taslim,  et  al.  (Ref  4)  calculated  the  flow  about  a  circular 
cylinder  and  a  109  thick  ellipse  in  pure  rotation  and  combined  rotation  and 
translation.  Their  stream  function-vortieity  formulation  was  based  on  the 
inertial  velocities  which  were  convected  with  the  body-axis  velocities  in  the 
body  fixed  reference  frame.  Their  calculation  used  an  explicit  method  based  on 
computing  one  velocity  component  at  a  given  point  by  integrating  the 
contributions  of  the  vorticity  in  the  exterior  flow,  on  the  body  surface,  and 
inside  the  body.  The  other  component  was  obtained  through  continuity. 
Pressure  gradients  at  the  body  surface  were  obtained  from  the  tangential 
component  of  the  momentum  equation.  Two  point  differences  were  used  in 
time,  upwind  differences  were  used  for  eonvective  terms,  and  diffusion  terms 
were  calculated  from  central  differences. 

The  approach  taken  in  this  work  is  similar  to  that  used  by  Galloway 
(Ref  5)  in  simulating  the  free  flight  of  a  flat  plate.  Galloway  used  the 
contravariant  velocity  formulation  and  the  same  numerical  method  as  Hodge. 
This  technique  was  used  to  calculate  the  trajectory  of  a  flat  plate  in  free  flight. 
The  free  flight  trajectory  was  a  marked  contrast  to  the  other  works  above, 
which  all  used  a  prescribed  trajectory  for  the  body.  Galloway  derived  the 
contravariant  velocity  relationships  by  starting  with  the  equations  written  for 
an  inertial  control  volume.  He  then  expressed  the  inertial  quantities  in  terms 
of  the  time  varying  transformation  from  the  inertial  to  the  body  frame,  and 
applied  the  chain  rule  to  obtain  the  derivatives  in  the  body  frame.  This 
resulted  in  equations  that  are  identical  to  the  present  method,  but  did  not 
provide  any  insight  as  to  the  physical  meaning  of  the  different  terms.  The 
Poisson  equation  for  pressure  was  used,  and  free  stream  Dirichlet  boundary 
conditions  were  employed  to  simulate  the  far  field. 

The  present  work  differs  from  that  of  Galloway  in  several  respects.  The 
governing  differential  equations  for  the  rotating  coordinate  system  are  derived 
from  first  principles,  rather  than  starting  with  the  differential  equations  for 
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the  unaccelerated  case. 


The  present  numerical  method  differs  from  Galloway’s  by  using  an 
arbitrary  grid  about  a  thick  body,  and  in  its  treatment  of  the  viscous  and 
pressure  terms.  A  generalised  transformation  on  an  "0"  type  grid  is  used 
instead  of  the  orthogonal  aHa  type  grid  that  Galloway  used  to  represent  the 
flat  plate.  The  pressure  and  viscous  calculations  are  based  upon  a  four  point 
aXa  configuration  central  difference  operator.  The  exterior  boundary  condition 
is  also  treated  differently  so  as  to  attempt  to  allow  the  disturbances  in  the 
wake  region  to  exit  the  computational  domain.  The  final  difference  between 
the  present  work  and  Galloway's  is  the  use  of  cell  centered  pressures,  rather 
than  collocating  the  numerical  values  of  the  pressure  at  the  same  grid 
locations  as  the  velocities.  Oborin's  method  for  solving  the  pressures  is  used 
throughout  the  domain,  instead  of  just  at  the  wall  as  was  done  by  Hodge, 
Hegna,  and  Galloway. 

The  principal  objective  of  this  work  was  to  demonstrate  the  use  of  the 
contravariant  velocity  formulation  for  an  arbitrary  body  and  arbitrary  motion. 
The  cell  centered  pressures  were  selected  to  avoid  the  difficulties  experienced 
by  the  previous  investigators  in  calculating  the  pressures  at  the  wall.  Chorin’s 
method  was  selected  over  the  Poisson  equation  for  its  simplicity  in  calculating 
the  pressures  at  the  cell  centers  from  velocities  known  at  the  corners  of  the 
cell.  The  differencing  scheme  for  the  viscous  terms  was  an  outgrowth  of  the 
differencing  used  to  calculate  the  pressure  gradients,  and  was  motivated  by  the 
desire  to  calculate  the  momentum  equations  in  the  strong  conservation  form. 
Subsequent  results  showed  that  the  nonconservation  form  of  the  convective 
terms  gave  better  results,  and  so  the  strong  conservation  form  was  used  for 
only  a  minority  of  the  cases.  Finally,  the  inviscid  boundary  conditions  were 
selected  in  an  attempt  to  model  the  wake  region  more  accurately  than  the  free 
stream  boundary  conditions  used  by  Galloway. 
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CONTINUUM  EQUATIONS 


The  momentum  equations  used  for  this  investigation  were  the 
contravanant,  or  "mixed*  form  of  the  incompressible  Navier-Stokes  equations. 
The  traditional  approach  to  computing  the  flow  about  a  rotating  body  is  to 
express  all  of  the  equations  in  terms  of  the  velocities  written  in  terms  of  the 
rotating,  body  fixed  reference  frame.  This  leads  to  the  introduction  of  a 
number  of  additional  terms  in  the  momentum  equations  because  Newton's 
Second  Law  is  only  valid  for  velocities  and  time  derivatives  taken  with  respect 
to  the  inertial  reference  frame.  These  additional  terms  are  often  represented 
as  fictitious  body  forces  (Ref  6),  so  that  the  conservation  of  momentum  can  be 
expressed  in  the  same  form  as  for  the  non-accelerating  case.  Even  though 
represented  as  forces,  these  terms  can  be  handled  as  part  of  the  momentum 
flux  and  are  treated  as  such  in  the  present  numerical  study  (c.f.  below). 

A  purely  inertial  description  of  the  momentum  equations  is  not 
convenient  because  the  position  of  the  body  is  constantly  changing  in  the 
inertial  frame.  This  would  preclude  the  use  of  a  fixed  grid  in  a  numerical 
calculation.  The  contravariant  velocity  approach  uses  inertial  velocities  which 
are  resolved  in  the  unit  vectors  of  the  body  fixed  coordinate  system.  This 
affords  a  considerable  simplification  of  the  equations,  as  they  are  very  similar 
to  the  momentum  equations  for  the  non-accelerating  formulation.  The 
fietictious  forces  of  the  traditional  approach  do  not  appear  explicitly  because 
they  are  implicitly  contained  within  the  inertial  velocity  components.  Galloway 
derived  the  equations  by  expressing  the  inertial  velocities  in  terms  of  the  time 
varying  transformation  to  the  body-axis  orientation,  and  substituted  them  in 
the  conventional  Navier-Stokes  equations  written  for  the  inertial  frame.  The 
approach  of  the  present  investigation  was  to  express  the  momentum  using 
inertial  velocities  and  describe  the  momentum  transport  in  terms  of  body-axis 
velocities.  The  resulting  equations  are  equivalent  to  Galloway's,  and  provide 


considerable  insight  as  to  the  physical  meaning  of  the  acceleration/rotationai 
terms  and  the  eigenvalues  of  the  momentum  equations. 

Kinematics 

This  section  defines  the  notation  used  to  describe  the  position  and 
velocities  of  points  on  the  body  and  in  the  fluid.  It  is  important  to  distinguish 
the  difference  between  the  rotating  body  fixed  reference  frame,  and  the 
unaccelerated  inertial  reference  frame  when  differentiating  vector  quantities 
with  respect  to  time.  The  notation  and  kinematic  definitions  are  taken  from 
Likins  (Ref  2).  The  notation  indicates  the  relative  relationships  between  points 
and  the  notation  for  time  derivatives  of  vector  quantities  explicitily  indicates 
relative  relationships  between  different  coordinate  frames.  Superscripts  to  the 
left  of  the  variable  name  indicate  that  the  time  derivative  was  computed  as 
seen  by  an  observer  in  that  reference  frame.  Subscripts  to  the  right  of  the 
variable  name  are  used  to  indicate  relative  relationships. 

For  example,  let  0,  P,  and  Q  represent  points,  with  0  being  fixed  in 
inertial  space,  Q  fixed  to  the  body,  and  P  located  anywhere  in  the  fluid,  as 
shown  in  Figure  1.  The  position  of  P  relative  to  Q  is  denoted  as  rP/q.  This 
position  vector  may  be  resolved  into  components  which  are  parallel  to  the 
instantaneous  unit  vectors  of  any  frame  of  reference  without  regard  to 
whether  the  unit  vectors  of  that  frame  are  rotating  with  respect  to  another 
frame  .  The  velocity  of  P  relative  to  Q  relative  to  the  inertial  reference  frame 
is  the  time  derivative  of  r P/q: 

‘drpfq  *  1^p/q 

<it  (  1) 

The  velocity  of  P  relative  to  Q  relative  to  the  rotating  body  reference  frame 


+  <L{  *Up/o|Yf  (  Up/o-*-RY-  Uq/o)  +  Xf  (  Vp/0“RX-lVq/0)|} 

dr} 
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a?  p  j  a*  drj 
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dr)  p  J  3£  dr)  (  48) 
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drj  p  J  3£  dr)  (  49) 

Notice  that  these  equations  may  be  written  in  the  form 

afu)  +  am  +  ak)  =  aih>  +  aii>  >  m 

at  df  drj  dr)  (  50) 

in  which  the  quantities  {u},  {f},  {g},  {h},  {i>,  and  {k>  are  column  vectors 

representing  each  of  the  terms  in  {},  in  turn,  for  both  equations  above. 

These  equations  have  substantially  the  same  form  as  those  written  for 
a  control  volume  fixed  in  the  inertial  coordinate  system.  The  difference  are 
the  {k}  vector,  which  results  from  differentiating  the  time  derivative  in  the 
body  frame,  and  in  the  mix  of  velocities  in  the  convective  terms.  For  the 
inertially  fixed  control  volume  ‘VP/q  is  a  constant,  and  therefore  may  be 

taken  to  be  xero  since  it  will  affect  none  of  the  terms  in  the  differential 


be  represented  as 


j  ac  =  a_(Y^  c)  -  a_(Yf  c) 
ax  dt  dr)  (  45) 


j  ac  =  -  a_(Xt,  c)  +  ajx*  c) 

ay  at  dr)  {  46) 

When  these  relations  are  substituted  into  J  times  u  times  the  Laplacian,  one 
obtains  the  result 

J  3_  (i/aC)  +  J  3_  (i/3C) 

ax  ax  ay  ay 

-  a_{  v  [aac  -  <sac  |>  -*■  a_{  v  [-/sac  +  73c  |> 

df  J  d£  dr)  drj  J  3£  dr)  (  47) 

with: 

a  -  (X,)2  +  (Yij)a 
fi  -  X*X„  +  Y*Y, 

7  -  (X*)a  +  (Yf)a 
J  »  IX^Y,  -  X,Y*| 

Notiee  that  no  cross  derivatives  appear  explicitly. 

The  final  step  in  the  derivation  of  the  continuum  equations  is  to  apply 
these  results  to  the  momentum  equations  after  multiplying  them  by  J  (which 
is  fixed  with  respeet  to  time)  and  dividing  by  p: 


3_{J‘upjo}  -*■  3_{  ‘up/oJY,  ('up/o-t-RY-'u^^o)  -  (1vp/o-"RX-tvq^o)|} 

at  at 


(  41) 


8*Vp  ft\  +  flJIUp/o1Vp;o-»-lVp;o(RY-1UQ/o)|  +  flJlVp/oa-,-lVp;o(-RX-*VQ;o)l 

at  ax  ay 

=  -  1  3fi  +  a3  *Vp /n  +  a2  Vp/n)  -  R  'up/o 

p3Y  ax2  3Y2  (  42) 

III  the  equations  above,  wb/iXbVp/q  has  been  moved  to  the  right  hand  side  as 
the  source  terms  R*vP/0  and  R‘up/0.  These  equations  are  identical  to  the 
equations  derived  by  Galloway. 


The  next  step  in  the  development  of  the  continuum  equations  is  to 
apply  a  general  transformation  from  the  physical  domain  to  the  computational 
domain.  £  represents  the  ordinate  and  rj  the  abscissa  in  the  computational 
domain.  The  continuity  equation  is  easily  transformed  using  the  chain  rule 
and  the  inverse  relationships,  Y,j/J=£x#  XTj/J=^Y)  Y^/J =tjx,  Xf/J=TjY.  as: 


d?.«PlQ  Ytf  -  dl\lv>/n  Yf  -  a*Vp/n  X,  +  3*Vp/n  =  0 
3£  drj  d$  dr)  (  44) 

X^,  Xq,  Yf,  Y jj  indicate  the  partial  derivatives  of  X  and  Y  with  respect  to  £ 
and  rj. 

In  order  to  put  the  momentum  equations  in  conservation  form,  the 
chain  rule  and  inverse  relations  above  are  first  applied.  As  shown  by  Viviand 
(Ref  8),  the  Jacobian,  J,  times  the  derivative  of  some  quantity,  C,  can  be  can 


rXY  3  UldlVyfn  + 

ax  ay 


Tn  -  M(aiVp/n  +  aWp^n) 

3Y  az 

TXZ  *  M(3*Up/n  +  aiVp/n) 

az  ax  (  38) 

For  incompressible  flow  in  which  the  viscosity  is  assumed  constant,  the 
right  hand  side  of  the  momentum  equation  becomes: 

RHSx  *  #1  -an  +  uia^  *up|n  +  a3  ‘up^n  +  aiSj^oJi  dv 

cv  ax  ax2  ay2  az2 

rhsy  *  #[  — an  +  ttfa2  ‘Tp/n  +  a2  *vp ^  +  a2  *vp^n)i  dv 
cv  ay  ax2  ay2  az2 

RHSj  *  ffi\  -dv  +  [Aid?  ‘wp^  +  a2  *wp^»  +  a2  1wp^d)[  dV 

07  az  ax2  ay2  az2  (  39) 

The  integrands  on  both  sides  of  the  equations  are  equated  for  each  component 
in  turn,  as  with  the  continuity  equation,  to  obtain  the  differential  equations. 

For  two  dimensional,  incompressible  flow,  the  flow  will  be  assumed  to 
lie  in  the  X-Y  plane.  Therefore,  P  «  Q  *  'wP/0  =  0.  The  angle  if  will  be 
used  to  represent  the  orientation  of  the  body  axes  relative  to  the  inertial  axes, 
as  shown  in  Figure  2,  and  its  rate  of  change  is  R.  Under  these  conditions  and 
the  continuity  and  momentum  equations  become: 

a*up|n  +  a*vp^n  ■  o 

ax  ay  (  40) 


ajsEjQ  +  aj‘up/o2+lup/o(RY~l,iq/o)l  -*■  3J‘up/olvp/o‘*,lup/o(-RX-lVQ/o)I 

at  ax  ay 
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M  (  ^(p'up/oj-^piQ^p/o-R'yp/o) 

cv  at 

+  V’jp'up/ot1  Vp^o-Wb/iXrp/Q-1  Vq/o)I  }  dV  =  RHSx 
M  ( ajp^p/ohpfR'up/o-P'vp/o) 

cv  at 

^p  vp/orVp/o-Wb/iXrp/q-'VQ/o)!  }  dV  =  RHSy 
M  (  ^_(p1'WP/o)',’P(PlvP/0_QluP/o) 

cv  at 

+  VV^P/oK'V p/o-^b/iXrp/g-'V q/0)]  }  dV  a  RHSj  (  35) 

The  forces  consist  of  pressure  and  viscous  forces.  The  effect  of  gravity 
on  the  fluid  will  be  neglected.  These  forces  are  typically  represented  as 
(Ref  7): 


with 


RHSx  ~  I\  -§£  +  (a  gy  +  a  t yy  +  a  t«m  dv 


cv 

ax 

ax 

ay 

az 

RHSy 

*  S\  -ao  + 

fa  Tw  ♦  a  try 

+  a  Tvsii  dv 

cv 

ay 

ax 

ay 

az 

RHSj 

*  m  -ai  + 

fa  tvi  +  a  Tvi 

-  «»  a  &x ))  dv 

cv 

az 

ax 

ay 

az 

°x  -  -2  fi  7.‘vP/0 

2/l^Ue^ 

3 

ax 

cry  =  -2  jx  7.lVP/0 

+  2udlyptn 

3 

ay 

<7j  a  —2  ft  7*‘Vp/0 

+  2/ia*Wp/n 

3 

az 

(  36) 


(  37) 


becomes: 


M I  !8(p'Vp,o)  *  x  (p'vP,0)|  dv 
ov  at 

^  p1  V  p/ol(l  V p/o- Wb/i^rp/q-1  V q/oHI  dA  »  RHS 
°s  (  33) 

This  equation  can  be  derived  from  the  conventional  fictions  body  force 

formulation  as  given  in  Shames  (Ref  6).  This  vector  equation  can  be 

expressed  in  terms  of  its  components: 


M\  &.(p‘ttp/o)+p(Q1vp/o-R,vp/o)l  dv 

07  at 

+  $  pV/oIPV p/O'^b/iXfp/Q^V (jyo)vi|  dA  *  RHSx 
cs 

JJF\  3_(plvp/o)-*-p(R1ap/o-P1Vp/o)l  dV 

w  at 

+  <$>  P^P/olPV p/O-^/i^rp/Q-1  V qyoHi  dA  ■  RHSy 
cs 

JR  9_(plvP/o)-*-p{P1vp/o-Q‘up/o)|  dV 
cv  at 

+  <$>  p'wp/ott’V p/o-«b/iXrP/Q-1  V q/o)t»I  dA  =  RHSj 
05  (  34) 

where  ‘up/0,  lvP/0,  ‘wP/0  are  the  inertial  velocity  components  in  the  X,  Y,  and 

Z  directions;  P,  Q,  R  are  angular  rates  about  the  X,  Y,  and  Z  axes;  and  RHSX> 

RHSy,  RHSi  are  force  components  applied  in  the  X,  Y,  and  Z  directions. 

When  the  Divergence  Theorem  is  applied  to  each  component,  in  turn,  the 

following  result  is  obtained: 
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acceleration  of  the  body-axis  reference  center,  Q,  does  not  affect  Vp/o>  The 
other  two  terms  represent  the  net  momentum  flux  out  of  the  control  volume 

If  one  describes  the  inertial  velocity  vector  field  using  the  time  varying 
position  vectors  of  the  body-axis  system,  the  integrals  for  region  III  and 
region  II  may  again  be  combined  by  recognising  that  point  P  is  being  used  as 
a  dummy  variable  and  evaluating  the  integrands  at  the  same  position  in  the 
body  axis  coordinate  system.  Once  again,  the  difference  divided  by  At  is 
represented  as  a  partial  derivative,  because  the  position  in  the  body-axis 
system  is  held  constant.  The  rotation  of  the  body  frame  with  respect  to  the 
inertial  frame  is  still  important,  however.  The  momentum  flux  through  the 
surface  of  the  control  volume  is  handled  as  with  any  other  extensive  property. 
Again,  there  is  no  flux  for  particles  in  the  flow  which  appear  to  be  moving 
tangential  to  the  boundary  in  the  body-axis  system. 

Applying  the  above  relationships,  as  with  the  continuity  equation, 
results  in: 


g  !<L(P1vp/o)  dV  -*■  S'  Ply p/o(by p/q^)  dA 
cv  at  cs 

-  #  f,  dA  -  g  fb  dV  -  RHS 

cs  cv  (  32) 

fa  is  the  net  force  per  unit  area  acting  on  the  surface  of  the  control  volume 

and  ft,  is  the  net  body  force  per  unit  volume.  RHS  simply  refers  to  the  right 

hand  side  (force  terms)  of  this  momentum  equation.  The  time  derivative  of 

the  momentum  within  the  control  volume  taken  with  respect  to  the  inertial 

frame  (which  has  unit  vectors  in  the  inertial  frame),  can  be  expressed  in  terms 

of  a  derivative  taken  relative  to  the  body  frame  (which  has  unit  vectors  in  the 

body  frame).  The  body-axis  velocity  may  also  be  expressed  in  terms  of 

inertial  velocities.  After  applying  these  changes,  the  momentum  equation 


differential  element  m  turn. 

The  impulse  applied  to  the  system  between  t  and  t+At  is: 

I  *  1  PlVP/0dV  -  Jjf  p'Vp/odV 

in  n 

+  1  ^Vp^dV  -  g  plVp/0dV 
IV  I  (27) 

The  impulse  is  the  integral  from  t  to  t+At,  taken  with  respect  to  the  inertial 

frame,  of  the  forees  applied  to  the  fixed  mass  system: 


I  -  /F«fc 


(  ») 


Differentiating  the  impulse  with  respect  to  time  vith  respect  to  the  inertisl 
frame  recovers  the  force  acting  on  the  system: 


(  29) 


To  first  order,  the  derivative  above  may  be  calculated  ar. 


)  -  Mft)  I 

At  | 


(  30) 


Applying  the  above  relationships  to  the  integral  expression  for  the  impulse 


results  in: 


F  -  {  IJVV^odV- JV Vp,0dV|/At  .  |jyV,,(0<tV|Mt 


-  I  JVV„odV|/At  (t 

1  (  31) 

The  first  bracketed  expression  represents  the  time  rate  of  change  of  the 
momentum  within  the  control  volume,  differentiated  with  respect  to  the 
inertial  frame.  The  distinction  that  the  momentum  is  differentiated  with 
respect  to  the  inertial  frame  accounts  for  the  changing  orientation  of  the  unit 
vectors  of  the  body  frame  with  respect  to  the  inertial  frame.  The  translational 


form  of  the  continuity  equation  is  also  the  same  as  for  the  non-rotating  case: 


V-‘VP/0  -  0  (21) 

Conservation  of  Linear  Momentum 

The  derivation  of  the  expression  for  the  conservation  of  linear 
momentum  for  an  accelerated  control  volume  proceeds  in  a  parallel  fashion  to 
that  for  the  the  conservation  of  mass.  Consider  again  the  fixed  mass  system 
which  occupies  regions  I  and  II  at  time  t  and  regions  III  and  IV  at  time  t+At 
The  force  acting  on  a  particle  at  a  given  time  is  equal  to  the  mass  times  its 
acceleration  relative  to  the  inertial  frame.  Thus  the  force  acting  on  a  particle 
containing  a  fixed  amount  of  mass,  dm,  and  located  at  point  P  is: 

dF  -  dm  Mf  rp/o  *  dm  M  *VP/0 

dt  dt  (  22) 

The  linear  momentum  of  the  particle  at  point  P  is  defined  as  the  product  of  its 
mass  and  its  inertial  velocity: 

dM  -  dm'Vp/o  (  23) 

The  force  acting  on  the  particle  is  the  time  rate  of  change  of  the  momentum 
of  the  particle,  because  the  constant  mass  can  be  brought  inside  of  the 
differentiation.  The  momentum  of  the  system  at  time  t  is  obtained  by 
integrating  the  momentum  of  the  particles  througout  regions  I  and  IL 

jBf  p  'Vp/o  dV  +  jj*  p  ‘Vp/0  dV 
I  II  (24) 

The  momentum  at  time  t+At  is  given  by: 

Jp'V P/0dV  -  Jp'Vp, 0dV 
m  IV  (  26) 

Point  P  is  used  as  a  dummy  variable  in  these  integrals  to  indicate  each 


Integrating  the  efflux  around  the  entire  eontrol  surface  and  putting  it 
into  the  expression  for  the  conservation  of  mass  for  the  fixed  mass  system 
yields  the  integral  statement  of  the  conservation  of  mass  for  the  accelerated 
eontrol  volume: 

Jif  d  p  dV  »  -$  p*V p/q*n  dA 

07  at  08  (16) 

Or,  expressing  the  velocity  in  terms  of  inertial  quantities: 

#  3  p  dV  *  p{iVp/oHtfb/iXrpfQ“,VQfo)*n  dA 

07  dt  (  17) 

Applying  Gauss's  Divergence  Theorem  yields: 

Jf  d  p  dV»  -jjf  Vip('V p/o-Wb/iXrp/Q-1  V q/o)1  dV 

^tat  07  (  18) 

This  equation  must  hold  for  all  control  volumes  with  a  fixed  volume,  even  as 

the  control  volume  is  shrunk  to  a  differential  sise,  therefore,  the  integrands  on 

both  sides  of  the  equation  must  be  the  same,  which  leads  to  the  differential 

form  of  the  conservation  of  mass: 

3p  »  -  V  •  (pOVp/0  -  Wb/i  X  rp/q  -  ‘Vq/o)) 
at  (  19) 

The  divergence  of  'Vq/o  -  0  since  it  is  constant  throughout  the  entire  control 

volume  and  only  varies  with  time.  The  divergence  of  «b/,  X  tP[q  is  sero, 

because  ub/i  is  also  spacially  constant  throughout  the  domain.  Thus  (he 

divergence  of  the  mass  flux  is  the  same,  whether  expressed  in  terms  of  the 

body-axis  velocity  or  the  inertial  velocity.,  and 

3p  *  -V-(p‘Vp/o) 

dt  (  20) 

which  is  the  same  as  the  classic  result  derived  using  a  control  volume  fixed  in 
inertial  space.  For  incompressible  flow,  the  density  is  constant,  and  the  final 


■JT 


as  At  goes  to  zero,  yields: 

Inn  £  loft+At)  -  of  til  dV  -  £  dp  dV 

At-X)  ov  At  cv  dt  (  13) 

where  the  partial  derivative  has  been  used  to  represent  the  instantaneous  time 

rate  of  change  at  a  given  point  fixed  in  the  body-axis  coordinate  system.  This 

is  appropriate  because  the  density  at  a  given  point  in  the  body  fixed  control 

volume  can  be  written  as  a  function  of  its  position  in  the  body-axis  system 

and  time,  ie,  p(rS/q,t).  It  is  rS/Q  which  is  being  held  fixed  for  this  partial 

derivative. 

The  net  mass  flux  out  of  the  control  volume  can  be  obtained  by 
integrating  the  efflux  through  each  element  of  the  control  surface.  The 
position  of  a  particle  in  the  system  can  be  written  using  its  position  relative 
to  the  point  on  the  control  surface  by: 

rP/0  *  *>/S  +  Ts/q  +  Tq/o  (  14) 

The  inertial  velocity  of  the  particle  then  becomes: 

‘Vp/O  =  'Vp/s  +  ‘Vs/q  +  ‘Vqyo 

lVp/o  *  bVP/s  >  «i/iXrp/S  +  bVS/q  +  Wh/iXrs/q  +  '^q/o  (  15) 

The  efflux  from  the  control  volume  at  point  S  consists  of  the  particles  that  are 
crossing  point  S  at  time  t.  The  efflux  though  the  differential  element  of  the 
control  surface  located  at  point  S  is  pbVp/s*n  dA,  where  n  is  the  outward 
normal  unit  vector  to  the  surface  at  S.  This  expression  is  the  result  of  the 
simple  observation  that  a  particle  located  on  the  boundary  at  time  t  which 
does  not  cross  the  boundary  must  be  stationary  or  moving  tangent  to  the 
boundary,  as  observed  from  the  body  frame.  Note  that  rp/s  *  0  at  time  t  for 
those  particles  on  the  control  volume  boundary,  and  ‘vs/,  -  0  by  the 
definition  of  the  body  fixed  control  volume.  Thus,  ‘Vp/s  -  “Vp*. 
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body-axis  system,  ie, 


"vs/,  -  0  t  «l 

and  there  is  a  one-to-one  mapping  of  ail  points  on  the  boundary  of  region  II 
to  the  boundary  of  region  III,  regions  II  and  III  may  be  regarded  as  defining 
a  body  fixed  control  volume.  Note  that  this  also  results  in  a  one-to-one 
mapping  of  all  of  the  points  within  the  tvo  regions  as  well  as  the  points  on 
the  boundaries. 

The  mass  in  the  system  at  time  t  is  given  by: 

g  p  dV  +  JpdV 

I  II  (9) 

and  the  mass  at  time  t+At  is  given  by: 

g  P  dV  +  g  p  dV 

III  IV  (  10) 

Since  the  system  has  the  same  mass  at  both  times,  there  can  be  no  difference 
between  the  two  expressions: 

#pdV  -  g  p  dV  +  jr  p  dV  -  g  p  AM  =  0 

III  II  IV  I  (  11) 

The  difference  between  the  first  two  integrals  above  indicates  the  total  change 

in  the  amount  of  mass  contained  within  the  control  volume  between  time  t  and 

t+At.  The  second  difference  represents  the  net  mass  flux  out  of  the  control 

volume  over  the  time  interval.  Because  of  the  one-to-one  relationship  of  the 

points  in  regions  II  and  III,  the  first  two  integrals  may  be  combined: 

J  p  dV  -  g  p  dV  -  g  |p(t-At)  -  p(t)|  dV 
IH  II  ov  (  12) 

provided  that  the  difference  in  density  at  the  two  times  is  evaluated  at  the 

same  point  in  the  body  fixed  axis  system.  Dividing  by  At  and  taking  the  limit 


relative  to  the  inertial  frame,  will  be  termed  the  "inertial  velocity"  of  a  particle 
in  the  flow.  This  may  be  regarded  as  the  "absolute"  velocity  of  a  particle  in 
the  fluid.  bYp/qJ  the  time  rate  of  change  of  the  position  of  the  particle  relative 
to  the  body-axis  reference  point  Q  as  seen  by  an  observer  rotating  with  the 
body  frame,  will  be  called  the  "body-axis  velocity." 

The  body-axis  velocity  is  calculated  from  the  inertial  velocity  by: 

bVp/q  »  ‘Vp/o  -  Utyi  X  rp/q  -  lVq/o  (8) 

This  quantity  is  seen  to  vary  with  the  translational  velocity  of  point  Q 
relative  to  0,  with  the  angular  rate  of  the  body  frame  relative  to  the  inertial 
frame,  and  with  the  position  of  point  P  in  the  body-axis  system.  The  inertial 
velocity,  ‘VP/0  is  utterly  independent  of  any  motion  of  Q  or  rate  of  rotation  of 
the  body  frame. 

The  only  quantities  whieh  will  be  referred  to  as  forces  in  this  paper  are 
those  forces  for  which  Newton's  second  law,  F  «  m  a,  holds  in  an  inertial 
reference  frame  Thus,  the  sum  of  the  forces  acting  on  a  particle  of  fixed 
mass  m  at  point  P  is 

F  =  m  rP/0  »  m  |d  Vp/o 

dt*  dt  (  7) 

Conservation  of  Mass 

The  derivation  of  the  relationship  defining  the  conservation  of  mass  for 
an  accelerating  control  volume  begins  by  examining  a  fixed  mass  system  of 
particles,  as  shown  in  Figure  2.  The  fixed  mass  system  at  time  t  consists  of 
the  particles  in  region  I  and  region  II  (regions  I  and  II  are  mutually  exclusive, 
as  are  regions  HI  and  IV).  At  a  later  time,  t+At,  the  fixed  mass  system 
occupies  regions  III  and  IV.  Let  point  S  represent  any  arbitrary  point  on  the 
boundary  of  region  IL  If,  at  t+At  point  S  maintains  its  position  in  the 
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is  expressed  as: 


dt  (2) 

The  rate  of  rotation  of  the  unit  ret  tors  which  define  the  orientation  of  the 
body  frame  relative  to  the  inertial  frame  is  denoted  by  and  the  velocities 
relative  to  the  inertial  and  body  frames  are  related  by: 

‘Vpjq  *  Mfp/q  *  +  wb/i  X  rP/q  *  bVP/q  -►  ttfc/i  X  TP/q 

dt  dt  (  4) 

Note  that  the  unit  vectors  which  define  the  orientation  of  a  reference  frame 
are  not  bound  to  any  point  Thus  the  only  time  rate  of  change  of  the  unit 
vectors  is  their  angular  rate,  u.  Translational  velocities  are  handled  by 
explicitly  noting  the  points  at  the  head  and  tail  of  the  position  vectors.  Thus, 
the  inertial  velocity  of  P  relative  to  O  is  related  to  the  velocity  of  P  relative 
to  Q  by: 

‘VP/0  -  ‘VP/Q  +  *Vqy0  (  5) 

For  the  rest  of  the  paper,  the  points  P,  Q,  and  O  will  have  a  special 
significance.  Point  0  will  be  taken  to  be  a  point  which  has  no  acceleration 
when  viewed  from  an  inertial  reference  frame.  Thus,  there  exists  some 
Galilean  transformation  for  which  the  point  0  may  be  considered  to  be  "fixed." 
Point  Q  will  refer  to  the  reference  center  of  the  body  fixed  axis  system. 
Positions  in  the  body-axis  system  will  be  measured  relative  to  point  Q.  Point 
P  will  be  used  to  designate  the  location  of  a  particular  particle,  having  a  fixed 
mass,  in  the  fluid.  rP/0  is  the  position  of  point  P  measured  from  the 
unaccelerated  point  0,  and  rP/q  is  the  position  of  the  particle  relative  to  the 
accelerating  body-axis  reference  point. 

‘VP/0,  the  time  rate  of  change  of  the  position  of  the  particle  relative  to 
the  unacceierated  reference  point  as  viewed  by  an  observer  who  is  not  rotating 
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equation.  With  ‘VP/q  -  0,  ‘Vp/0  ■  lVP/q  =  bVP/q.  The  terms  in  parentheses 
then  become  equal  to  ‘Vp/Q,  {k}  goes  to  sero,  and  ‘Vp/Q  can  also  be  substituted 
for  ‘Vp/0,  thus  eliminating  the  mixed  use  of  velocities  from  different  reference 
framse.  The  advantage  of  describing  the  flow  in  terms  of  the  inertial 
velocities  results  from  the  simplicity  of  the  terms  arising  from  the  time 
derivatives.  The  source  term,  {k},  is  half  the  magnitude  of  the  typical  Coriolis 
acceleration.  The  centripetal  acceleration  and  the  angular  acceleration  do  not 
appear  explicity.  This  is  an  advantage  numerically,  as  these  terms  depend 
upon  the  magnitude  of  rP/q,  which  approaches  infinity  at  the  exterior  of  the 
computational  domain.  The  acceleration  of  the  body  is  also  not  required.  The 
disadvantage  in  using  the  inertial  velocity  is  a  loss  in  the  intuitive 
interpretation  of  the  flow  about  the  body.  The  body-axis  velocity  is  easily 
calculated  from  the  inertial  velocity,  however,  and  so  this  disadvantage  is  not 
great. 

Convective  Term  Eigenvalues 

The  numerical  approximation  to  the  fluxes  will  depend  on  eigenvalues 
of  the  isolated  convective  terms.  These  are  most  conveniently  found  by 
considering  the  equations  in  linearised  form.  By  applying  the  chain  rule, 
Equation  31  becomes: 

aim  +  Faim  +  qaim  =  afh>  +  am  +  m 

at  df  dr)  df  drj  (  51) 

where: 

{U}  *  J‘Vp/0 

f  -am 

3{U> 
g  -aiii 


21 


The  matrices  F  and  G  are 


F  = 


Yjj(bup/Q+‘up/o)  -  XqbvP/q  -Xjj‘up/o  1 

J  J  J  I 


Y,lvP/0 

J 


Xnb«p/Q  -  Xj^vriq+'vpto) 
J  J 


(  52) 


G  - 


-X^(bttP/g'*’,«P/o)  +  X^Vp/g  2$4'uP/0 

J  J  J 


-Y$bUp/q  +  X$(bvp/g+iVp/o) 
J  J 


(  53) 


with 

bUp/g  »  1Up/0-'-RY-1Ugyo 

bYP/g  -  ‘Vpyo-RX-Vgyo 

If  continaity  is  exactly  satisfied,  then  the  nonconseryatiye  form  is  equally 
valid,  and  the  matrices  F  and  G  become: 


1 


Ygb«p/q  -  Xqbyp/q 

»  9 


«p/q  +  X4  vp/q  0 
J  J 


(  54) 


0  -XtSg-XtV.,, 

(  55) 

These  matrices  are  also  the  same  as  would  be  obtained  by 
quasi-linearising  the  equations,  assuming  that  the  body-axis  velocity  is  locally 
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constant  This  assumption  is  better  tor  regions  of  the  flow  which  are  some 
distance  from  the  body,  so  that  the  body's  disturbances  are  small  compared 
with  the  magnitude  of  the  "freestream*  flow.  Since  the  body-axis  velocity  is 
the  negative  of  the  inertial  velocity  of  the  grid  plus  the  inertial  velocity  of 
the  flow,  the  inertial  velocities  may  be  viewed  as  perturbations  to  the  grid 
velocity. 

The  eigenvalues  of  F  and  G  (assuming  continuity  is  satisfied)  are: 

\t  *  Yj^up/q  -  Xqbyp/q 
J  J 

Xq  »  ~Y$bUp/q  +  X^Vp/q 

J  J  (56) 

The  eigenvalues  are  seen  to  depend  upon  the  body-axis  velocity,  not  the 
inertial  velocity.  This  is  a  direct  consequence  of  the  use  of  the  body-axis 
velocity  to  describe  the  transport  of  momentum  through  the  control  volume. 
The  sign  of  these  eigenvalues  indicates  the  flow  of  information  due  to  these 
terms.  The  numerical  algorithm  depends  upon  the  sign  of  these  eigenvalues 
to  determine  in  which  direction  the  upwind  finite  differences  will  be  made. 

This  result  indicates  an  important  rule  for  understanding  the  role  of  the 
body  axis  and  inertial  velocities.  Whenever  one  is  using  a  velocity  in  a 
context  which  regards  the  convective  transport  of  any  quantity,  the  appropriate 
velocity  to  use  is  the  body-axis  velocity.  If  the  context  regards  the  linear 
momentum  contained  within  the  flow,  inertial  velocity  is  the  appropriate  one. 
For  example,  if  one  is  calculating  the  stream  function  by  integrating  along  a 
line  in  the  computational  grid,  the  body-axis  velocity  should  be  used.  This  is 
because  the  stream  function  represents  the  mass  flux  accross  the  line  used  in 
the  integration.  Only  in  this  manner  will  the  stream  function  have  a  constant 
value  at  the  body  surface. 

Near  the  body,  the  body-axis  velocities  are  dominated  by  the  inertial 


velocity,  since  at  the  body  surface  the  inertial  velocity  is  identical  to  the 
velocity  of  the  grid  and  the  body-axis  velocity  is  zero  (no  slip  boundary 
condition).  Ideally,  continuity  should  1m  satisfied  at  all  times,  but  this  may 
not  be  true  numerically  as  the  solution  is  iterated  to  convergence  at  a  given 
time  step.  If  this  is  the  case,  then  the  additional  terms  of  the  conservative 
form  may  be  a  source  of  error  in  the  momentum  equations.  Therefore,  it  is 
useful  to  examine  the  continuity  terms  as  they  appear  in  the  the  momentum 
equations.  The  differences  between  the  conservation  and  nonconservation 
forms  of  F  and  G  are: 

XquP/o  -Xjj'up^o 
J  J 

Xq'vP/o  ~Xql*p/o 

J  J  J  (58) 


G’  - 


(  59) 


If  F’  and  G’  are  combined  along  with  their  associated  velocity  derivatives,  the 
result  is  {u}(V*‘VP/0).  The  eigenvalues  for  F’  and  G’  are: 


Aft  3  XqUp/o  “  Xq‘vp/n 
J  J 

Ap*2  3  0 

Ao't  *  “X^p/o  *  X$‘vp/o 

J  J 

A<«  3  0  (  60) 

The  associated  eigenvectors  for  the  non-zero  eigenvalues  are  (U)  for  both 
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matrices.  If  continuity  is  not  satisfied  due  to  incomplete  numerical 
convergence,  the  upwind  differencing  of  these  terms  must  be  considered, 
especially  if  conservative  differencing  is  used.  For  any  area  of  the  flow  in 
which  the  body-axis  velocities  are  decelerated,  the  eigenvalues  of  F  and  F’  or 
of  G  and  G’  may  be  of  opposite  sign.  This  happens  close  to  the  body,  in 
separated  regions,  and  near  the  stagnation  streamline.  These  error  terms  will 
then  be  downwind  differenced  and  may  be  unstable. 


The  continuum  equations  are  approximated  numerically  with  finite 
difference  equations.  The  finite  difference  equations  are  then  solved  using  a 
block  successive  over  relaxation  (SOR)  method.  Each  column  vector  in 
Equation  31  is  differenced  so  as  to  model  the  physical  dynamics  which  it 
represents.  Second  order  special  differences  are  used  throughout  For  the 
remainder  of  this  paper,  U  and  V  will  be  understood  to  be  ‘up/0  and  lvP/0  (the 
inertial  velocities). 

A  curvilinear,  body  fitted  grid  is  used  to  define  the  mapping  to  the 
computational  domain.  This  is  shown  in  Figure  2.  The  computational  domain 
is  a  rectangular  area,  in  which  the  grid  lines  are  straight,  orthogonal,  and 
spaced  one  unit  apart  (A£*l,  Arj-l).  The  geometry  of  the  body  and  grid,  and 
the  velocities  are  defined  at  the  mesh  node  points,  and  indexed  by  i  in  the  £ 

direction  and  by  j  in  the  tj  direction.  The  metric  quantities  Xf,  Yf,  X^f  Y^, 

and  J  are  defined  at  the  center  of  the  quadrilateral  ceils  bounded  by  the  mesh 
points.  In  addition  to  the  metrics,  the  pressure  (P),  viscosity,  and  velocity 
gradients  are  also  defined  at  the  ceil  centers.  The  cell  centers  are  defined  in 
the  transformed  plane  rather  than  the  physical  plane,  and  the  geometrical 
coordinates  of  the  centers  are  not  explicity  calculated.  The  center  points  are 
indexed  similarly  to  the  mesh  points,  ie.,  the  point  (i+l/2,j-*-l/2)  lies  within  the 
cell  bounded  by  (i,j),  (i+lj),  (i,  j+1),  (i-*-l,j+l).  These  indices  will  be  indicated 
using  subscripts,  and  the  time  level  of  the  variables  will  be  indicated  using 
superscripts.  The  index  n  will  indicate  the  present  time,  and  n-1  will  indicate 
one  time  step  (At)  in  the  past 
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The  velocity  components  ire  differenced  in  time  using  %  first  order 
backward  difference. 


m.  -  ( mi  -  mra  )/At  -  ( « -  «■*  )j/At 

dt  (  61) 

This  results  in  in  implicit  method  with  no  stability  limitations.  The  time  step, 
At,  is  therefore,  only  constrained  by  the  accuracy  required  to  adequately  model 
the  unsteady  flow  dynamics. 

The  convective  terms  are  differenced  using  three  point  one  sided 
differences  which  are  differenced  in  the  upwind  direction.  The  upwind  direction 
is  determined  by  the  sign  of  the  eigenvalue,  as  calculated  by  Equation  36.  The 
upwind  differences  are,  for  positive  eigenvalues: 

m  =  3/2  {t)l  -  2  {f} f.tiJ  -•  1/2  {fKU, 

as 

31a1  «  3/2  {g}*  -  2  {g>“H  *  1/2  {gft  _2 

drj  (  62) 

and  for  negative  eigenvalues: 

m  -  -3/2  mi,  -  2  {f>f+lil  -  1/2  miu, 

as 

m  -  -3/2  {g)s:(  -  2  {g)“i+1  -  1/2  {g}“+s 

3  r)  {  63) 

One  sided  differences  were  selected  for  these  terms  because  they  correctly 
modelled  the  flow  of  information  represented  by  the  fluid,  and  because  they 
have  dissipative  truncation  error  characteristics  that  enhance  the  stability  of 
the  calculation.  Metric  terms  appearing  within  the  flux  vectors  were  calculated 
using  conventional  central  differences. 

A  nonlinear  consistency  analysis  was  performed  for  these  terms 
differenced  in  both  conservative  and  nonconservative  linearised  form.  The 
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uonconservative  form  is  obtained  by  evaluating  the  diagonal  F  and  G  matrices 
at  the  i,j  location  for  all  terms.  This  is  in  contrast  to  differencing  them  at 
each  point  in  the  mesh  molecule,  as  is  done  implicitly  when  differencing  the 
{f}  and  {g}  vectors  in  the  conservative  manner  above.  For  simplicity,  the 
analysis  assumed  a  uniform  rectangular  grid  and  differencing  such  that  AX  => 
A£  and  AY  =  Arj.  For  the  conservative  differencing,  this  results  in: 

3_(1«p/obup/q)  3_(luP/obvp/q) 
dr) 

-  V/qUf  +  ^u^Vtj)  +  bvP/Qu1f 

-  [(2uf-*-R)uff  +  l/3tu-*hup/Q)u^i(A4-)2 

*  l^Yp/qU^  +  l/3uV^](A77)a 

*  OMfJAvft  (  W) 

4.(‘vP/ob«p/q)  +  i_(lVp/obvp/q) 

%  rj 

*  bUp/qVf  ♦  v(uf+v,)  +  bVp/QV?j 

-  ((Tfiiff  ♦  (uf-*-R)vff  +  l/3vu^f 

+  l/Shup/qVfffKAf)3 

-  I(2v1,-R)viw  ♦  l/^v-fhvp/qJv^KArj)2 

*  (  «) 
The  unadorned  velocities  and  their  derivatives  are  inertial  velocities  at  the  i,j 


grid  point.  Using  nonconservative  differencing  for  the  same  case  results  in: 


bup/o3Jup/o  +  bvp/o£Jup/o 
%  V 

-  buP/qUf  +  brP^uv 

-  |l/3buP/qufffl{A£)a  -  [l/3bvP/gu^{(^U)a 

+  OKAtrjAvr\  (  66) 

bup/o5_‘vp/o  +  hyp/o^'vp/o 

W  V 

=  V /qVf  +  brp/gTij 

-  (l/3bUp/QVf^l(A|)a  -  [l/3bVp;Q)Yjpp}](A7})2 

+  OKA^^A^l  (  67) 

Both  methods  are  second  order  accurate.  When  compared  with 
nonconservatire  differencing,  conserrative  differencing  is  subject  to  a  diffusion 
truncation  error  whose  sign  depends  on  the  sign  of  the  rotation  and  the  first 
order  velocity  derivatives,  and  thus  may  be  stabilising  or  destabilising  in 
different  parts  of  the  flow..  It  also  includes  the  velocity  divergence  terms,  as 
expected.  These  latter  terms  are  identically  sero  in  the  ideal  case. 
Numerically,  however,  these  terms  may  not  be  identically  sero  due  to  many 
reasons.  These  include  not  solving  directly  for  all  of  the  points  in  the  domain, 
but  updating  them  in  a  point  wise  or  row  wise  manner;  errors  in  the  pressure 
calculation;  and  incomplete  convergence  of  the  equations  in  general. 
Nonconservative  differencing  may  be  expected  to  be  more  robust  under  these 
conditions  because  of  the  fewer  truncation  error  terms. 


The  pressure  gradient  and  viscous  terms  were  differenced  using  a  four 
point  central  difference  scheme.  This  seheme  results  from  a  continuous  linear 
interpolation  operator,  which  was  based  upon  the  ideas  of  Walitt  (Ref  5).  The 
linear  interpolator  for  some  quantitiy,  X,  with  i£*£i-l,  is  given  by: 

i  X-ti+WXH-’})  *  Xf^-iKH-rjJ 
-  X^i-l-fX^H)  -  X.V.tf-ifo— j)  (  68) 

This  function  reduces  to  linear  interpolation  at  the  boundaries,  and  at  the 

center,  it  becomes  the  average  of  the  corner  points.  If  the  derivatives  with 

respect  to  *  and  rj  are  also  evaluated  at  the  center,  they  too  become  the 

simple  average  of  the  central  differences  at  the  edges: 

Xf  (H+1/2.TH+1/2)  »  (  -x&  -  X?+liJ  -  Xt%  *  X,V.  )/4 

X?|  tf-i-l/2, rj-j-l/2)  -  (  -X&  -  Xf+I|i  -  X“i+i  -  X“+tiJ+i  )/4  (  69) 

The  operator  is  second  order  accurate,  as  seen  from  the  modified  equations  for 
the  *  and  rj  derivatives: 

Xf  (f *  Xf  «-  Xfff(Af)V24 

-  Xf„(A»j)7*  *  0((Afr,  (At,)4| 

X,  -  X,  ■*  X„,IAt,|724 

*  Xff,(Af )7*  ■.  0[(Af  )4.  (At))4]  (  71) 

This  central  difference  operator  is  used  repeatedly  to  calculate  the  (h)  and  {i} 
vectors  at  the  cell  centers  by  calculating  the  metrics  and  velocity  gradient 
terms  from  the  surrounding  mesh  nodes,  as  in 
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ufi+1/2,,+1/2  =  (  +  «i+l,j  -  <+,  +  )/* 

■*•1/2, 1+1/2  *  (  “  U?+t.,  +  <+i  +  »?+t.j+l  )/4 

V^i+t/2,,+1/2  »  (  ”<i  +  Vf+l,,  -  V“j+l  +  Yf+,,+1  )/4 

v^+1/2,1+1/2  -  (  -  vf+t|J  +  <j+t  +  v?^  )/4  (  72) 

in  order  to  form  the  {h}  and  {i}  vectors.  Once  the  (h)  and  (i)  vectors  are 
calculated  at  the  cell  centers,  the  four  point  central  difference  operator  is  again 
applied  to  calculate  their  derivatives  at  the  mesh  node  point  of  interest  The 
pressure  gradient  at  the  node  point  is  also  calculated  from  the  cell  centered 
pressures  in  the  same  manner. 

Discrete  Momentum  Equations 

When  the  finite  difference  approximations  above  are  substituted  into 
the  momentum  equations  (28-29)  the  convective  portion  becomes  (for  positive 
eigenvalues): 


{J/At  *  3/2HY„„  -  Yt J  k<  *  (Xf,„  -  X,„)  kv”|)  <  - 

-  2KY?,-,,  ‘of-,.,  -  x;Wj  kv*,„)  of.,, 

/«■#■  K  n  «aifl  K  ft  \  Ml 


kiyi..,,,  y, 

—  Yn  bv: 

,1  Aqi-t,i  Ti 

•*  (Yfij.,  y, 

*  Xfi,,.,  kT 

-  l/2KY^.,., 

y-,,  -  x;,. 

-  ( -  x?,,-, 

y,.,  -  xf,,r 

-  (J/At)  u“-? 

+  RHSX 

-2,1 

-2) 


(J/At  *  3/JKY,,.,  -  Y(,,,j  -  (Xf„  -  X*,)  kv“]}  vj,  *  R  < 

-  «YU  V-,.,  -  XV,., 

k<,.,  *  Xf,.,.,  k*»H)  r“p,| 


'  >*ll  *  Ml  •  *  Ml  W  Ml  " 

-  l/2((Y^-2tJ  bu?_2t)  -  X{j,_2il  bvf.2il)  vf.2iJ 

(••n  k  n  mm  k  n  v  n 

-  I  fi 


b,.n 


V! 


with 


bu?i  *  ut\j  +  R  Yii  “  ‘u<3/o»“i 

H  =  <>  -  RX,i  -  Vq/o^ 


The  rotational  source  terms  have  been  brought  back  to  the  left  hand  side 
above,  in  preparation  for  the  block  SOR  calculation.  In  contrast  to  the  {h}  and 
{i}  vectors,  the  metric  terms  in  the  {f}  and  {g}  vectors  are  calculated  using 
conventional  central  differences  about  the  node  points: 


Xfij  -  (Xwj  -  Xi.J/2 
X^ij  =  (Xy*.  -  X,iH)/2 
Y{„  -  (Y..,,  -  Y..J/2 
Y,„  -  (Y.j.,  -  Y,  rl)/2 


(  H) 


Calculating  the  metrics  at  each  i,j  location  and  differencing  them  with  the 
fluxes  allows  an  intuitive  interpretation  of  these  terms.  The  centrally 
differenced  metrics  may  be  considered  to  represent  the  face  areas  of  a  cell 
which  is  centered  about  the  grid  line.  Multiplying  this  area  by  the  velocity  at 
the  line  is  an  approximation  to  the  total  momentum  flux  moving  through  that 
face.  These  fluxes  are  then  differenced  to  calculate  the  net  flux  moving 
through  the  cell.  This  is  analogous  to  a  finite  volume  approach. 

The  force  terms  become: 
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RHS*  *  l/2{  [P/p(Yf  -  Y,)|  i+1/2, j+1/2  ~  [P/p(Yjj  >  Yf)|f+l/2,_1/2 
-  [P/rfY,  -  Yf*W^  +  !p/M  Y*  +  } 

+  l/4{  ufi  [  -  »f+i/2, 1+1/2  -  *f-l/2, j+1/2  -  *f-l/2,|-l/2  “  ^f+1/2,1-1/2 
+  2(bf<.1/2ij+1/2  -  bf_1/2iJ>i/2  -  bi+1/2i).1/2  -*■  bi.1/2r,/2) 

“  gf+l/2,,-1/2  ~  gf*t/2,j-i/2  “  gf-l/2,i+l/2  “  gf-t/2, 1-1/2! 

*  uf+l,!+l  (*f+t/2.i+t/2  ~  2bf>1/2|+1/2  4-  gf+i/2,i+l/2l 

-*■  «f+l,,-t  (»f+l/2, 1-1/2  +  2bf^l/2  rl/2  +  gf+1/2,1-1/21 

+  uf-,,1+1  [*f-l/2, j+1/2  +  2b,“  t/2il+i/2  +  gi— 1/2,1+1/zl 

+  uf-,,,-1  (ftf.|/2,j-t/2  -  2bui/2  j+l/2  4-  gf-1/2, 1-1/21 

+  af+1,1  [»f+l/2, 1+1/2  +  *f+l/2.j-l/2  ~  gf+1/2, 1+1/2  ~  g  1*1/2, 1-I/2I 

+  **“,+1  [gl+1/2, t+1/2  +  gf-1/2, j+1/2  ~  lf+l/2, 1+1/2  ~  li-l/2, 1*1/2! 

+  ttf-l,i  (*f-l/2,j+t/2  *f-t/2,j-t/2  “  gf-1/2, j+1/2  “  gf-t/2,|-l/2l 

+  Ufri  [gf-1/2, 1-1/2  +  gf+i/2,|-t/2  “  *f-t/2, j-l/2  “  ^+1/2,1-1/21  } 

RHSy  =  l/2{  [P/p(Xtj  -  )|j+i/2, 1+1/2  ~  [P/p(Xij  +  Xf)|f-i/2, ,+1/2 

+  [P/p(X*  -  X,)|f.,/2(rl/2  4.  [P/p(X,  4.  X*)|f>,/2, r,/2  } 

+  l/4{  v“  [  -  &f»l/2,j+i/2  -  *f-l/2,j+l/2  ~  *f-l/2, 1-1/2  “  *f+l/2, 1-1/2 
*  2(b,“  1/2, ,+i/2  -  bf.,/2>1vi/2  ~  bf*i/2ir(/2  bf_,/7r(/2J 

-  gf+l/2,j-l/2  ~  gf+1/2, j-l/2  “  gf-1/2, 1+1/2  ~  gf-1/2, 1-1/2! 
vf-l,i+l  [*f-t/2,j+l/2  *  2bf_,/2  j+l/2  4-  gf-1/2,, +1/2! 

+  vf-t,i-l  W-l/2, 1-1/2  ~  2bf^t/2|J+,/2  4-  gf-t/2,,-i/2l 
+  vf+t,l+l  W+l/2, 1+1/2  ~  2bf+t/2j+,/2  4-  gf.,/2, 1+1/2! 

+  vf*l,|-l  (*f+l/2, j-l/2  +  2bf+1/2iJ_,/2  4-  gf+,/2,j_i/2l 
+  vf+l,j  (*f+l/2, j+1/2  +  ^f+1/2, i-t/2  “  gf+1/2, j+1/2  ~  gf+1/2, 1-I/2J 

vfi+l  [gf+l/2,j+l/2  gf-t/2, j+1/2  ~  *f+l/2, j+1/2  “  *f-l/2, j+l/2l 

+  Tf-l,i  !*f-!/2, 1+1/2  *f-l/2,j-t/2  “  gf-t/2, 1+1/2  ~  gf-1/2, J-I/2I 

+  Tfj-l  [gf-l/2, j-l/2  +  gf+l/2,j-t/2  ~  ^f- 1/2,1- 1/2  ”  *{‘+1/2, j-l/2l  }  (  75) 

where: 

a  =  i/a/J 
b  =  v(S/J 
g  =  vy/J 

Again,  the  quantities  at  i±l/2,  jdtl/2  are  defined  at  the  centers  of  the  cells 
suirounding  the  i,j  node  point. 
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The  nine  point  combination  of  velocities  representing  the  viscous 
contribution  is  still  second  order  accurate.  As  an  example,  consider  a  straight, 
uniform,  but  not  necessarily  orthogonal  grid  in  which  a,  y,  and  fS  are  constant. 
For  this  case,  the  viscous  terms  (RHSviiJ  become: 

RHSvucx  =  (-*  -  g)  <  +  (*  -  g)  Uui.j/2  +  (a  -  g)  uf_tiJ/2 
-*■  (~»  +  g)  <*i/2  *  {-a  +  g)  u“j_,/2 
+  (»  -  2b  +  g)  u^.,^,/4  +  (a  +  2b  +  g)  u?_,iH.,/4 
+  (»  +  2b  +  g)  u“.£j.1/4  (a  -  2b  +  g)  uf.,,r,/4  (  76) 

The  other  component  has  the  same  form.  When  Taylor  series  are  substituted 

into  the  equation  above,  the  modified  equation  is  obtained: 

RHSvucx  *  a  uff  -  2  b  uf,  +  g  uvv  +  a(A£)2  u^/12 

-  MA^)2  Uf^jj/3  +  (g(A£)2  +  a(A7j)2]  utfvv/4 

-  b(AV)*  uhl„/3  -  g(Arj)2  u^/12  *  0((A£)4,(A7j)4l  (  77) 

The  cross  derivative  term  has  appeared  implicitly  even  though  it  was  not 
expicitly  written  into  the  original  difference  equation. 

Next,  these  equations  are  solved  for  u“  and  v“.  This  requires  the 
solving  for  them  simultaneously,  since  the  rotation  terms  couple  the  two 
components  together.  The  equations  may  be  put  into  the  form  [A]{u}={B), 
where  (A|  is  a  two  by  two  matrix  and  {B}  is  a  column  vector  containing  all  of 
the  terms  except  for  the  ones  containing  u“  and  v“.  The  elements  of  [Aj  and 
B)  are  (for  positive  eigenvalues): 

Atl  =  Aw  -  J/At  -  3/2KY,  -  Yfa  bu*  +  (X*  -  X,),,  bv“] 

1/ 4|a,“.1/3j+1/2  >  if-l/h+t/l  •+■  lf»t/2,rl/2  +  *?-l/2, 1-1/2 
-  2(b“l/2  ~  ^f-l/2,1+1/2  “  ^f+l/2,1-1/2  ■+■  ^i-l/2, 1-I/2) 

I?*l/2,|-*l/2  gf-J/2,1+1/2  +  gf*l/2, 1-1/2  8?-l/2, 1-I/2] 

A,2  *  -  R 

A*  -  R  (  78) 
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CONCLUSIONS  AND  RECOMMENDATIONS 


The  objective  of  developing  a  method  to  calculate  the  laminar 
incompressible  Navier  Stokes  equations  about  arbitrary  bodies  undergoing 
arbitrary  motion  was  not  completely  achieved. 

The  derivation  of  the  Navier  Stokes  equations  for  an  accelerating 
control  volume  was  accomplished.  The  nonconservative  and  strong 
conservation  form  for  these  equations,  including  viscous  terms  in  conservation 
form,  was  derived.  The  use  of  inviscid  boundary  conditions  allowing 
disturbances  to  exit  the  domain  was  also  investigated,  but  not  conclusively 
demonstrated. 

The  numerical  solutions  have  been  characterized  by  slow  convergence 
and  difficulties  in  conserving  mass  throughout  the  domain.  There  is  some 
evidence  to  suggest  that  there  are  long  period  periodic  error  components 
which  are  not  being  eliminated  by  the  relaxation  method  used  to  solve  the 
equations. 

It  is  recommended  that  further  investigations  be  conducted  to  resolve 
these  problems.  Approaches  which  might  be  profitable  include: 

1.  Abandoning  the  cell  centered  pressures  and  artificial  compressiblity 
method  in  favor  of  the  more  traditional  Poisson  equation  for  pressures 
evaluated  at  the  nodal  points. 

2.  Using  a  different  method  to  solve  the  discrete  equations.  Possible 
candidates  would  include  using  a  conjugate  gradient  algorithm  to 
minimise  some  measure  of  the  residuals,  and  multigrid  methods. 

3.  Investigating  different  boundary  conditions,  particularly  for  the  cell 
centered  pressures. 
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The  pressure  contours  highlight  the  asymmetric  nature  of  the  solution. 
They  also  indicate  that  there  may  be  some  difficulty  in  the  forward  stagnation 
region.  The  change  in  slope  crossing  the  periodic  boundary  is  an  indication 
that  a  similar  error  to  that  seen  in  Stokes  first  problem  may  be  occuring. 

The  remaining  two  frames  show  the  global  domain.  Once  again,  the 
major  effects  of  the  solution  are  confined  to  the  immediate  vicinity  of  the 
body,  indicting  a  lack  of  convergence  in  the  flovfield. 
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free  stream  condition  for  this  case.  This  is  also  shown  by  the  streamlines, 
pressure  contours,  and  vorticity  contours. 

The  remaining  plots  show  the  same  information,  but  for  the  entire 
domain.  Some  pressure  influence  has  propagated  throughout  the  domain  in  a 
symmetric  fashion.  The  looping  shape  of  the  pressure  contours  suggests  that 
there  may  be  a  periodic  error  which  has  a  period  of  one  half  the  width  of  the 
computational  domain. 

The  third  test  case  is  the  same  circular  cylinder  which  has  been 
impulsively  accelerated  from  rest  to  a  constant  speed  at  a  Reynolds  number  of 
20.  .  At  the  time  of  the  results  shown  in  Figure  8,  the  cylinder  has  moved  a 
distance  of  0.4  diameters.  Ten  time  steps  have  been  used  with  20  SOR 
iterations  at  each  time  step.  The  inviseid  boundary  conditions  have  been  used 
at  the  far  field. 

The  inertial  velocity  vectors  of  Figure  8a  indicate  that  very  little 
acceleration  of  the  flow  is  taking  place  at  the  midpoints  of  the  cylinder  and 
the  wake  region  shows  little  deceleration  of  the  flow  except  in  the  immediate 
vicinity  of  the  cylinder.  The  flow  is  also  not  symmetric.  These  observations 
are  also  indicated  by  the  body-axis  velocity  vectors  of  the  next  plot.  This  plot 
also  indicates  that  the  flow  is  attached  all  the  way  to  the  aft  stagnation  point. 

The  third  view  shows  the  streamlines.  The  stream  function  was 
calculated  by  integrating 

*  ■  /  Pb «p/q  dy  -  /  pbvP/q  dx  (  93) 

along  the  radial  grid  lines  extending  from  the  surface  out  to  the  extenor 
boundary.  The  streamlines  are  lines  of  constant  4r.  The  trapexoidal  rule  was 
used  for  integration.  A  careful  comparison  of  the  body-axis  velocity  vectors 
and  the  stream  function  contours  shows  that  the  velocity  vectors  are  not 
parallel  to  the  contours,  particularly  in  the  region  above  the  cylinder.  This 
indicates  that  mass  is  not  being  conserved  in  this  region. 


program  is  having  difficulty  in  eliminating  errors  in  this  area.  During  the 
initial  stages  of  the  calculation,  the  magnitude  of  the  velocity  was  seen  to 
decrease  across  the  entire  domain  except  at  the  very  last  point.  Here  the 
velocity  was  set  equal  to  the  velocity  at  the  other  side  of  the  domain  because 
of  the  periodic  boundary  condition.  The  error  in  the  velocity  thus  had  a 
sawtooth  pattern.  As  the  calculation  was  iterated,  the  mean  value  of  the 
velocity  at  a  given  height  increased,  and  the  magnitude  of  the  sawtooth 
decreased.  The  pressure  contours  are  reflecting  this  error  in  the  velocity 
calculation. 

The  final  plot  from  Stoke’s  first  problem  shows  contours  of  constant 
vorticity.  As  expected,  they  are  parallel  to  the  surface.  Since  the  ideal 
pressure  gradient  is  lero,  boundary  layer  theory  indicates  that  the  velocity 
profile  next  to  the  wall  should  be  linear.  Thus  the  normal  derivative  of  the 
vorticity  should  be  zero.  The  wider  spacing  of  the  vorticity  contours  next  to 
the  wall  indicate  that  this  condition  is  being  approximately  satisfied. 

The  next  case  to  be  considered  is  a  circular  cylinder  in  pure  rotation. 
Figure  6  shows  the  grid  used  for  this  case.  It  consists  of  46  points  in  the 
circumferential  direction  and  41  points  in  the  radial  direction.  The  radial 
distribution  is  geometrically  stretched  so  as  to  maintain  a  constant  aspect  ratio 
for  each  cell.  The  outer  boundary  is  at  25  radii  The  rate  of  rotation  was 
ehosen  to  give  a  Reynolds  number  (based  on  the  diameter)  of  100  at  the 
cylinder  wall.  Dirichlet  free  stream  boundary  conditions  were  used  again  for 
this  case.  100  time  steps  have  elapsed,  with  five  SOR  iterations  being 
performed  at  each  time  step. 

The  inertial  velocity  vectors  are  shown  in  Figure  7a.  It  can  be  seen 
that  the  solution  has  hardly  progressed  from  the  region  in  the  immediate 
vicinity  of  the  cylinder.  This  is  also  shown  by  the  body-axis  velocities  shown 
in  the  next  plot  The  no-slip  boundary  condition  is  enforced  and  the  flow 
rapidly  transitions  to  the  apparent  solid-body  rotation  which  represents  the 


Two  time  steps  were  used,  with  100  SOR  iterations  at  each  time  step. 
A  single  iteration  consisted  of  relaxing  each  point  along  a  line  of  ^constant 
starting  from  the  wall  and  proceeding  outward  to  the  outer  boundary.  This 
sweep  was  applied  to  each  vertical  grid  '.ine,  moving  from  left  to  right  across 
the  entire  domain.  It  was  then  repeated  while  moving  the  line  sweeps  from 
right  to  left,  also  sweeping  from  the  wall  to  the  outside  along  each  line.  A 
relaxation  parameter  of  1  was  used.  The  elapsed  time  was  chosen  so  that  the 
top  of  the  domain  at  the  final  time  corresponds  to  T)  =  2,  where  17  is 
Schlichting's  nondimensional  height,  rj  =  Y/(2>/yt).  The  plate  is  moving  from 
left  to  right. 

The  inertial  velocity  vectors  shown  in  Figure  5a  indicate  the  unknowns 
for  which  the  equations  are  solved.  Plotted  for  comparison  is  the  exact 
solution.  It  can  be  seen  that  the  qualitative  agreement  iB  good.  The  flow  is 
parallel  to  the  plate  and  the  profiles  at  each  location  are  the  same.  There  is 
some  numerical  disagreement  This  may  be  due  to  the  small  number  of  large 
time  steps  and  the  attendant  lag  caused  by  the  truncation  error.  It  may  also 
be  due  to  a  lack  of  complete  convergence  in  the  calculation.  The  velocity 
profile  was  still  evolving  slowly  when  the  program  was  terminated  after  100 
iterations. 

The  next  view  in  Figure  5  shows  the  body-axis  velocity  vectors.  These 
show  the  boundary  layer  profile  and  illustrate  the  no-slip  boundary  condition 
at  the  wall.  The  streamlines  presented  in  the  next  view  are  also  consistent 
with  the  velocity  vectors,  indicating  parallel  flow. 

Figure  5d  presents  lines  of  equal  pressure.  The  spacing  between  the 
lines  corresponds  to  a  change  in  cP  of  1.6x10" 7 .  This  verifies  that  for  the  outer 
bounday  conditions  chosen,  the  pressures  are  essentially  constant.  Over  most 
of  the  domain,  the  pressure  gradient  is  also  zero.  However,  near  the  corners 
of  the  domain  the  large  relative  changes  in  the  pressure  indicate  that  the 
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COMPUTATIONAL  RESULTS 


The  numerical  experiments  which  have  been  performed  indicate  that  the 
method  is  promising,  but  not  yet  fully  developed.  Test  cases  which  have  been 
tried  include  Stoke's  first  problem,  a  circular  cylinder  in  pure  rotation,  a 
circular  cylinder  in  pure  translation,  an  elliptical  cylinder  in  pure  translation, 
and  an  elliptical  cylinder  in  pure  rotation.  The  results  from  the  first  three  test 
cases  will  be  presented  to  illustrate  some  of  the  problems  that  have  been 
encountered. 

Stokes  first  problem  consists  of  a  quiet  fluid  bounded  by  a  flat  plate 
which  is  infinite  in  extent  in  both  directions.  The  plate  is  instantaneously  set 
in  motion  at  a  constant  speed.  This  results  in  a  shear  layer  next  to  the  plate 
which  diffuses  upward  as  time  progresses.  The  exact  solution  for  this  case  is 
presented  in  Schlichting  (Ref  7).  The  pressure  should  be  zero  throughout  the 
domain,  and  the  flow  is  parallel  to  the  plate.  In  addition,  the  velocity  profiles 
at  each  location  at  a  given  time  are  identical.  This  case  was  selected  as  a 
benchmark  because  it  has  an  exact  solution  for  comparison.  It  is  sensitive  to 
the  pressure  calculation  and  the  flow  is  parallel  to  the  outer  boundary,  which 
provides  a  difficult  test  of  the  exterior  boundary  conditions. 

The  computational  grid  for  this  example  is  shown  in  Figure  4.  It  is  a 
uniform  rectangular  mesh  with  16  points  in  the  £  direction  and  41  points  in  the 
rj  direction.  Periodic  boundary  conditions  are  applied  to  the  ends  of  the 
domain  (£=1  and  £*16).  Since  this  is  a  simply  connected  physical  domain, 
rather  than  the  doubly  connected  physical  domain  surrounding  a  closed  body, 
the  periodic  boundary  conditions  can  be  regarded  as  simulating  an  inifinite 
series  of  domains  stretching  to  the  right  and  left  of  the  domain  un4er 
consideration  at  a  given  iteration.  For  the  results  shown  in  Figure  S,  free 
stream  Dirichlet  boundary  conditions  were  applied  to  the  top  of  the  domain. 
Thus  this  calculation  is  equivalent  to  non-steady  Couette  flow. 
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along  the  entire  boundary.  Note  that  this  is  similar  to  applying  Dirichlet 
boundary  conditions  at  phantom  grid  locations  just  outside  of  the  boundary. 
This  had  a  negligible  effect  for  translating  cases,  since  the  eigenvalues  on  the 
downstream  side  of  the  grid  resulted  in  differencing  into  the  domain. 

Pressures  were  assumed  to  be  xero  outside  of  the  domain.  This  was 
modelled  after  the  common  experimental  observation  that  the  static  pressure 
rapidly  returns  to  the  free  stream  value  as  one  moves  away  from  the  body. 
Since  the  outer  boundary  was  typically  more  than  ten  chords  away  from  the 
body  surface,  this  was  judged  to  be  a  reasonable  assumption.  It  also  has  the 
effect  of  fixing  the  value  of  the  pressure,  since  the  absolute  value  of  the 
pressure  is  arbitrary,  and  only  its  gradient  affects  the  fluid  dynamics.  The 
effect  of  the  assumed  pressures  outside  the  domain  on  the  pressures  inside  is 
felt  in  an  indirect  way.  The  outside  pressures  only  affect  the  velocities  on  the 
boundary.  These  velocities  are  then  used  to  compute  the  pressures  at  the 
centers  of  the  cell  just  inside  the  boundary. 

No  means  were  explicitly  employed  to  ensure  that  the  net  mass  flux 
through  the  entire  boundary  was  conserved.  It  was  assumed  that  the  pressure 
calculation  would  enforce  the  conservation  of  mass  throughout  the  domain. 
This  is  the  principal  deficiency  in  these  boundary  conditions.  It  is  especially 
apparent  for  symmetric  cases  involving  flow  parallel  to  the  boundary,  such  as 
the  case  of  the  circular  cylinder  in  spinning  motion.  Here  there  is  no 
constraint  to  ensure  that  the  velocities  at  the  boundary  are  directed  in  a 
strictly  tangential  direction. 
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terms.  Viscous  terms  were  neglected  at  the  boundary  since  mathematically 
they  are  elliptic  in  nature  and  require  values  outside  of  the  domain.  This  is  in 
keeping  with  the  observed  fact  that  the  Euler  equations  are  a  good 
approximation  to  the  flow  outside  of  the  boundary  layer.  If  the  fluid  was 
flowing  into  the  domain,  the  upwind  difference  formula  required  inertial 
velocity  values  outside.  These  were  assumed  to  be  iero,  corresponding  to 
undisturbed  fluid.  If  fluid  was  exiting  the  domain,  the  upwind  differencing 
only  required  values  inside  and  at  the  boundary  itself,  and  no  boundary 
condition  was  required  for  these  terms.  The  determination  as  to  whether  fluid 
was  entering  or  leaving  the  domain  was  made  by  examining  the  sign  of  the 
eigenvalues  of  the  convective  terms. 

Two  variations  were  tried  for  the  outer  boundary  conditions.  In  the 
first,  the  eigenvalues  were  computed  in  the  same  manner  as  in  the  rest  of  the 
domain,  that  is,  based  upon  the  body-axis  velocities  at  the  boundary  points. 
The  second  method  was  similar  to  the  first,  except  that  the  inertial  velocities 
were  ignored.  That  is,  bVp/q  was  assumed  to  be  approximately  ~v<J/o  for  the 
purpose  of  determining  the  direction  of  the  upwind  differencing.  This  is 
generally  a  good  approximation,  since  the  exterior  boundary  should  be  far 
enough  from  the  body  that  the  inertial  velocities  will  be  small.  For  most 
cases,  there  was  a  negligible  difference  between  the  two  methods.  The  second 
method  was  somewhat  better  in  the  event  that  the  exterior  flow  iB  essentially 
parallel  to  the  entire  outer  boundary.  An  example  of  such  a  situation  is  the 
flow  about  a  c  icular  cylinder  in  pure  rotation  using  a  polar  grid.  If  the  flow 
had  a  slight  outward  component,  the  undisturbed  condition  existing  outside  the 
domain  had  no  influence  at  all  on  the  flow  inside  the  domain  when  the  first 
method  was  used.  There  was  nothing  to  fix  the  value  of  the  flow  at  the 
boundary.  With  the  second  method,  the  algorithm  was  biased  slightly  by 
differencing  outward  when  the  eigenvalue  was  exactly  iero.  Thus,  for  cases 
such  as  the  rotating  cylinder,  the  differences  were  in  the  outward  direction 


Three  types  of  boundary  conditions  are  required  for  this  problem.  The 
computational  domain  is  bounded  at  the  bottom  by  the  solid  body,  at  the  top 
by  the  fluid  outside  the  computational  domain,  and  on  each  side  by  the 
opposite  side  of  the  domain. 

At  the  surface  of  the  body,  the  no  slip  condition  is  applied.  This 
requires  that 

kvp„  -0  (91) 

When  applied  to  the  inertial  velocities,  this  requires  that  lVP/0  »  lVq/0,  where 
point  Q  is  taken  to  be  the  point  on  the  body  at  which  the  boundary  condition 
is  being  applied.  In  other  words,  the  velocity  of  the  fluid  must  follow  the 
velocity  of  the  wall. 

The  periodic  boundary  conditon  at  the  sides  of  the  domain  results  from 
making  a  branch  cut  from  the  ouside  of  the  physical  domain  to  the  body  in 
order  to  map  the  doubly  connected  physical  domain  to  a  singly  connected 
computational  domain.  This  boundary  condition  was  enforced  by  taking  points 
from  the  opposite  side  of  the  domain  when  the  difference  formula  would  have 
required  points  outside.  The  points  at  £  *  1  and  £  =  iMAx  were  at  the  same 
location,  and  represented  opposite  sides  of  the  branch  cut.  For  example, 
values  at  l-l  in  the  finite  difference  stencil  evaluated  at  i=l  were  obtained 
from  iMAX~l-  1°  addition,  the  values  at  i  =  iMAX  were  set  equal  to  those  at  i 
»  1  as  soon  as  they  were  computed. 

The  boundary  conditions  at  the  outer  boundary  were  designed  to  allow 
undisturbed  fluid  (lVP/o=0)  to  enter  the  domain  where  the  edge  of  the  domain 
is  advancing,  and  to  allow  disturbances,  particularly  the  momentum  defect  in 
the  wake,  to  exit  the  domain  where  the  edge  is  retreating. 

This  was  handled  through  the  upwind  differencing  of  the  convective 


boundary  conditions. 

An  alternate  pressure  calculation  to  the  formula  above  was  also 
investigated.  For  an  inertial  control  volume,  consider  the  case  in  which  the 
velocity  divergence  is  locally  constant  over  an  entire  ceil  in  the  computational 
grid.  For  this  case  pJV-V  represents  the  net  mass  flux  out  of  the  cell  The 
equation  above  can  then  be  interpreted  as  relating  the  pressure  within  a  cell 
to  the  net  production  of  mass  within  that  cell.  The  net  mass  flux  may  also  be 
calculated  by  integrating  around  the  boundary  of  the  ceil.  This  suggests: 


AP  -  -2  0  {J/  [At(a  -  y)|}  A* 


(  88) 


A*  -  $  pbuP/q  dy  -  /  pbvP/Q  dx  (  89) 

If  trapesoidal  integration  is  used  on  each  side  of  the  cell,  the  discrete 
approximation  is  obtained: 

AP?+i/2,j+i/a  ■  -0  p  /  fAt(aj+1/2il+t/2  7i*-i/2.j«-i/2)l 

»  (V.,.,.,-VI)(X,.,.-X,1.1) 

*  (b<,., ..XX.^.-X-J  |  (  90) 

This  method  converged  faster  than  the  first.  It  was  also  possible  io  use  the 
conservation  form  for  the  momentum  equations  for  this  method,  whereas  the 
original  one  led  to  instablities  unless  the  nonconservation  form  was  used. 


calculated  (as  in  Gauss-Seidel  iteration),  and  so  storage  at  both  iteration 
levels  is  not  required.  The  calculations  above  are  iterated  at  the  same  time 
step  until  sufficiently  converged. 


The  velocity  calculations  above  depend  upon  the  use  of  the  pressures  to 
enforce  the  conservation  of  mass  while  iterating  to  convergence.  The  artificial 
compressiblity  approach  was  selected  for  the  pressure  calculation.  In  this 
method,  the  change  in  pressure  between  iterations  is  proportional  to  the 
velocity  divergence  (Ref  6).  This  approach  has  been  shown  by  Hodge  (Ref  7) 
to  be  equivalent  to  the  SOR  solution  for  the  Poisson  equation  for  pressure  for 
inertial  control  volumes.  Hodge  also  presents  the  optimum  acceleration 
parameter,  ♦,  for  this  calculation  as: 


*  -  Q  J*  /  [At  [a  +  7)| 


(  86) 


AP  -  -  p*V-‘VP/0  (  86) 

When  applied  to  the  pressure  calculation,  this  results  in: 

AP?+i/2,j+j/2  =  -0  Ji*l/2,jfl/2  /  [At  (Cti«.j/2,j«-l/2  7i+l/2,j+l/2)l 

[f'^q>+>/3,i*J/2",''^'fi-*-l/2,j-*l/2)  ^i.j  ■+■  (Y£i+ij2,j«.i/2— 1/2) 

( Y^,+l/2ii+t/2- Y fi+i/2,j+i^)  u“i+l  -  ( Y  1}i+i/2,)+i/2-*- Y ^i+i/5, ,4.1/2)  uf^t 

+  1+1/2)  *“j  (X^i+,/2,j+i/2“X^j4.t/2,,^|/j)  vf+l(i 

(X^i*i/2,|+l/2-X^i+1/2,J+|/2)  V®J+i  -  (X^,/2,j+|/2+X^i+l/2,j*|/j)  V?*,^,  ] 

(  87) 

at  the  cell  center.  One  advantage  of  this  approach  is  that  the  pressure 
calculation  itself  requires  no  boundary  conditions,  since  the  pressure  is  only 
calculated  at  the  cell  centers  which  are  wholly  within  the  control  volume. 
Pressure  boundary  conditions  are  still  required  for  the  momentum  equations, 
but  these  may  be  be  easily  implemented  for  both  Dirichlet  and  Von  Neuman 
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B,  *  l/4(af.|/2l|4.|/2  +  *?-t/2,j-l/2  ~  8f-l/2,j+l/2  -  S»a-t/2,|-t/2l  vf-l,l 

-*•  1  /4j-*f-i/2,j-|/2  -  *f+l/2,j-t/2  +  Si— l/2,i— 1/2  +  8?+l/2,j-l/2]  vi!rl 
+  l/2|Yj)i+2jbu?«.2tj  -  Xtji+2ilhvf>2J  vf+2ii 

-*■  l/21Y^,ii+2bu,“*2  -  Xtjiil+2bv“+2]  v“j+2 
+  -  X^1+libvf+lj) 

+  l/4(*f+l/2.ifl/2  +  *i+l/2.H/2  “  S 1+1/2, i+l/2  ~  Sf+1/2,1-  1/2)1  Vf+t.j 
+  [2(-YfiiMh<j+1  + 

+  1/ 4(-ftf+i/2ll+i/2  -  *f-l/2,j+|/2  +  Sf+1/2, 1+1/2  +  S 1-1/2, 1+1/2)!  vy+l 
+  l/^(*i+l/2,i+l/2  —  *^i+l/2,j+l/2  "*■  Si+l/2, 1+I/2)  ^i+lj+l 
+  l/4|s(L|/2,j+|/2  +  2bf_i/2(j+i/2  +  gf-1/2, 1+I/2]  vf-l,j+l 
+  l/4(*f+i/2,i-i/2  -*•  2bi^i/3ij_i/2  +  81+1/2,1-1/2!  ▼i+lj-l 
■*  l/4(*£-i/2,j-i/2  ~  2bf_,/2  j.1/2  +•  81-1/2,1-1/2]  vf-i,j-i 

+  l/2{  [P / p(X|j  -  X^)|ui/2J+l/2  “  |P/P(X,  +  Xf)|f_i/2,1+l/2 

+  (P/MXf  -  Xij)]“-i/2,i-1/2  +  IP/MX,  +  Xf)|f+i/2ij_i/2  } 

+  CJ/At)  v*£  (  82) 

The  equations  above  have  been  grouped  so  that  if  sufficient  storage  exists,  the 
metric  quantities  in  brackets  may  be  calculated  once  during  grid  generation  and 
stored  for  use  during  the  flowfield  calculation. 

Solving  for  u“  and  v“: 

u*ti  *  (Bj  A22  ~  Ai2  B2)  /  (Ah  A22  -  A«  A21) 

v*5!i  *  (B,  Au  “  A2j  Bj)  /  (Ah  A 22  -  At2  A2i)  (  83) 

The  block  successive  over  relaxation  algorithm  (SOR)  first  calculates  the 
velocities  at  a  point,  as  above.  It  then  over  relaxes  the  velocities  using  the 
relaxation  parameter,  Q,  as: 

=  0  <  +  (1-0X1-" 

</"  -  0  <  +  (1  -  0)  </-"  (  14) 

where  the  superscript  (s)  indicates  the  present  iteration  level,  and  (s-1) 
indicates  values  at  the  previous  iteration  level.  0  must  be  less  than  2  for 
stability.  The  computer  algorithm  replaces  the  values  at  i,j  as  they  are 


>  > 


For  negative  eigenvalues,  the  equations  become: 


A„  -  A„  -  J/il  -  3/2((Y,  -  Y(J„  k<  *  (X,  -  X,)„  ‘y“| 

*  t  n  n  n 


B,  = 


l/4(af U/2,,+1/2  ■+■  *T-l/2,|+l/2 
“  2(h?+i/2,i+t/2  “  bf-t/2,1+1/2  “  ^f+1/2,1-1/2 

A  n  n 

+  Si+l/2,1+1/2 


,a 

Si-l/2,j+l/2 


*?+l/2,i-l/2  *?-i/2,|-l/2 

®  '  *>f-l/2, 1-I/2) 

I?+t/2,j-l/2  +  81-1/2,1-1/2) 
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-  R 
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1/4W  -l/2,j+l/2  **•  *?-l/2,i-t/2  “  «?-l/2,+„?  ”  4?-l/2,H/*l  U?-l,i 

+  1/ 4|— afL1/ajj_t/2  -  nf+i/2,j-i/a  8?-t/2,i-i/2  "**  8?+i/2,j-i/2l  u?j-i 

+  l/2)Y,l+a>uf.2il  -  X^j+2ilhvf>2j)  ut“2|i 

l/2)Yjjl]+2bu“>2  -  X1jiij^2bv“j4.2)  u^a 

*  (-2(  Y jji+i,,buf+lil  -  Xfji«.tjbv£.,J 

+  l/4(a|t,/2i^.,/2  +  af^t/^j-j/a  -  gf*i/2,j+i/2  -  lf+i/2, 1-1/2)!  «f+i,j 

*  l2(-Y^iiWbu,Vi  +  X^)1+,bv“i+i) 

*  l/4(”*£.i/2,j+i/2  “  *?-l/2,j+l/2  +  Sf+1/2,1+1/2  **■  gf-l/2, 1+1/2))  «i\i+l 
+  1/ 4(a{^i/2i)«.l/2  -  2bf+i/2j+,/2  ■+■  gf+i/2,,+i/2)  u?+i,j+i 

*  l/4(af_i/2|j+i/2  ■+•  2bf.,/2il<ll/2  -*■  gf-t/2|i+i/2)  uf_i>J+, 

.  1  /iL&  a.  OL&  ^  vn  1  A 

ai+i/2(j-l/2j  ui«.|,]-t 
a  1  n 


l/4)aj  ♦1/2, 1-1/2  +  2bf+1/2rl/2  +  g,“*t/2, 1-1/2)  «f+ij- 

*  1/ 4)af_i/2ri/2  -  2b“_,y2  j_,^2  -*■  gf-i/2,i-t/2)  uf-1,1-1 

*  l/2{  |P/p(Yf  -  Y -  IP/PtY,,  -  Yf*W: 

*  |P/p(Y,  -  Y*)]f_l/2iJ_1/2  +  jP/p(Y,  +  Yt)|f.,/JtJ+,/2  } 


1/2 


!  JfW«fcu  - 

*  I/4W  -l/2,j+l/2  +  *?-l/2,i-l/2  ~  4?-t/2,]+l/2  ~  4?-l/2,j-l  /2JI  uf-l,j 

-  I  *  Xf, 

+  l/4(“» f-|/2,i-l/2  ~  *?+l/2,j-l/2  +  4?-t/2.i-t/2  Si+t/2, 1-1/2)!  U,“-i 

-  l/2tYijj_2ilbu£.2,j  ~  uf-2,j 

-  l/%sk<,  -  X1ji)J_,bv“1_2l  <_2 

+  1/ 4(-»f+i/2iJ+i/2  -  X?+l/2,]-t/2  +  4?+l/2,j+t/2  Si^t-l/2, j— 1/2]  ai>l,i 

+  1/ 4|»?*i/2(,+i/2  +  *?-l/2,j+l/2  ~  4?+l/2,j+t/2  ~  gf-1/2, 1+1/2]  «?j+l 
+  1/ 4(*f+i/2, 1+1/2  ~  2b?*,/2tj*,/2  +  41+1/2,1+1/2]  ^i+lj+1 

*  1/4W  -t/2,1+1/2  2b?_1/2iJ*l/2  -*■  sf-1/2,i+l/2l  ttf-l.j+1 

+  l/4(»f+1/2,j-i/2  ■+•  2bf+1/2  j_l/2  ■+•  41*1/2,1-1/2]  U?+1,1-1 
+  1/ 4l»f-l/2,|-l/2  “  2bf_t/S|_,/2  +  4i— l/2,j— I/2I  uf-l.|-l 

-  l/2{  [P/p(Yf  -  Y  tj)lf+l/2, j+l/2  “  IP/plY,  -  Yf)|».l(!.H„ 

+  [P/p(Y^  -  Yf)|jLi/2j-i/2  +  [P/p(Yi}  +  Yf)|f.1/2,J+i/2  } 

-  (J/At) 

liKY^V-u  -  X„-, 

h-  l/4(i; 

-1/2, 1+1/2  ■*"  *?-t/2,j-l/2  ~  4f-l/2,i+l/2  “  4?-1/2,i-1/2)1  V?-l,j 

-  [2(-Yfitr,yr,  +  X^j-jhyUj.,) 

-*■  If 4(-»f_i/2,j-l/2  “  ^1+ 1/2,1- 1/2  +  4?-l/2,|-l/2  +  41+1/2,1-1/2)!  v!j-l 

-  ~  ^^i-2,ibvf-2,il  vf-2,i 

-  l/2jYVii-2b<-2  ~  Xtjiij_2bv“j_2]  v“_2 

+  1/ 4(»i+l/2,i+l/2  +  *?+t/2,j-l/2  ~  4?+l/2,i+l/2  ”  4 1+1/2, 1-I/2]  vf+l.i 

+  1/ 4(_i+t/2, ,+1/2  ~  *?-l/2,j+|/2  +  41+1/2,1+1/2  4f-l/2, 1+1/2]  v?j+l 

+  l/4]»f+|/2,i+l/2  ~  2bf*,/2,,*t/2  -*■  4l+t/2tj+l/2l  Vf+l,i+l 

-  1/4W  -t/2,1+1/2  +  2b“_1/2, 1*1/2  +  4f-l/2, 1+I/2]  vf-l,i+l 

-*■  l/4|xf*l/2(1  -i/2  +  2b“*1/2, 1-1/2  41+1/2,1-1/2]  Tf*i,i-i 

-*■  1/ 4{Xi-i/2,j-i/2  “  2bf_i/2,i_l/2  +  4f-l/2, 1-1/2]  Vf-!,1-1 

-  l/2{  [P / p(X^  -  Xf)|f*l/a,i*,/2  -  [P/piX,  -  X^lf.,/2, ,+,/2 

+  lP/p(Xf  -  ^1))ll— l/2,j-l/2  +  [P/p(Xij  +  X^)|f*i/2,i-i/2  } 

-  mi)  <,! 


(  79) 
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Figure  2.  Control  Volume/Control  Surface  Relationships 


Global  View 


Figure  4  a.  Computational  Grid-  Stoke’s  First  Problem 


Streamlines 


CONTOUR  SPACING  IS  . 50000E-0I 


Figure  5c.  Stoke’s  First  Problem 
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Vorticity  Contours 
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Figure  5e.  Stake's  First  Problem 
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Inertial  Velocity  Vectors 


Figure  7a.  Circular  Cylinder-  Pure  Rotation 


Body- Axis  Velocity  Vectors 


7b.  Circular  Cylinder-  Pure  Rotation 


Pressure  Isobars 


Figure  81.  Circular  Cylinder-  Pure  Translation 
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Vorticity  Contours 


Figure  8e.  Circular  Cylinder-  Pure  Translation 
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Pressure  Isobars 
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Streamlines 


Figure  9c.  Circular  Cylinder-  Pure  Translation 


Vorticity  Contours 
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Streamlines 


Figure  7h.  Circular  Cylinder-  Pure  Rotation 
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Body— Axis  Velocity  Vectors 
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Inertial  Velocity  Vectors 


Figure  7f.  Circular  Cylinder-  Pure  Rotation 
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Streamlines 
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Vorticity  Contours 


Figure  8g.  Circular  Cylinder-  Pure  Translation 
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— '  -ne  purpose  of  this  project  was  to  develop  a  numerical  method 
capable  of  calculating  the  laminar  flow  about  a  two  dimensional 
accelerating  body  using  the  incompressible  davier-Jtokes  equations. 

This  problem  arises  in  the  calculation  of  the  dynamic  stall  of  an 
airfoil  and  in  the  calculation  of  pitching  and  heaving  stability 
derivatives  for  an  isolated  airfoil.  It  is  also  a  first  step  toward 
solving  the  three  dimensional  flow  about  a  complete  wing  undergoing 
arbitrary  motion  at  highly angles  of  attack,  including  departure,  spin, 
and  post-stall  maneuvers  .v 

The  contra variant  f'orm  of  the  momentum  and  continuity  equations 
are  derived  and  a  finite  difference  approximation  developed.  The 
c on tra variant  velocity  formulation  is  based  on  the  use  of  both  inertial 
and  body-axis  velocities  to  achieve  a  form  for  the  momentum  equations 
that  has  the  potential  for  improved  numerical  characteristics  for 
rotating  coordinate  frames. 

Numerical  calculations  to  date  have  not  been  satisfactory. 

Test  cases  presented  are  3toke's  first  problem,  a  circular  cylinder 
translating  at  a  Reynold's  number  of  20,  and  a  rotating  stationary 
circular  cylinder.  . 
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